NASA-CR- 194392 


Annual Report 


//a */G -tr#*' 


Studies on Equatorial Shock Formation 
During Plasmaspheric Refilling 

Grant #NAGW-2128 
P.I.: Dr. N. Singh 




During the grant period starting August 1, 1992, our major effort has been on 
examining the presence of equatorially trapped hot plasma on plasmaspheric 
refilling. We performed one-dimensional PIC simulations of cold plasmas 
expanding into a hot plasma, consisting of hot anisotropic ions and warm isotropic 
electrons, trapped in a region of minimum magnetic field. Simulations showed that 
the electric potential barrier built up by the anisotropy of the hot ion population 
facilitates in the formation of electrostatic shocks when the cold ion beams begin to 
come into contact with the hot plasma. The shock formation occurs even when the 
cold ion beams are highly supersonic with respect to the ion-acoustic speed. This 
findin g is interesting because equatorial shock formation during the early stage of 
plasmaspheric re fillin g has been debated over about two decades. In the past ion- 
ion instability has been invoked as the main mechanism for the coupling between the 
cold ion beams approaching the equator from the conjugate ionspheres. This 
coupling occurs when the beams are sufficiently slow; the beam velocity being less 
than three times the ion-acoustic speed. In the presence of hot plasma, the beams 
slow down by the potential barrier. The slowing down and the reflection process 
lead to the formation of the electrostatic shock even for highly supersonic ion 
beams. 


The mixing of hot and cold plasma was also studied. It was found that the cold 
electrons eventually short-circuit the potential barrier and the cold plasma punches 
throu gh the equatorially trapped hot plasma. This spatial mixing produces velocity 
distribution functions which are potentially unstable. The perpendicular velocity 
distribution of ions consists of a cold core and a hot ring; the latter starts at the ion 
energy corresponding to the magnitude of the initial potential barrier set up by the 
plasma. Such ring distributions are known to be unstable exciting ion-Bemstein 
modes, as well as lower hybrid waves. These waves can in turn heat the cold core. 
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The results from this study were reported in a Huntsville Workshop (October 
15-17, 1992) on plasma models, AGU 1992 Fall meeting [Singh, 1992] and a paper 
written for Geophysical Research letters [Singh, 1993]. An extended paper for JGR 
is under preparation. 

During this grant period we have also continued our investigation of the 
comparison of the large-scale hydrodynamic and kinetic models for plasmaspheric 
r efillin g comparison brings out the shortcomings and the strengths of the respective 
models [Singh, 1993]. This study also reconciles some of the previously reported 
results from these models. 
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Abstract . Effects of equatorially trapped hot plasma on 
the highly supersonic cold-plasma flow occurring during 
early stage plasmaspheric refilling are studied by means 
of numerical simulations. It is shown that the equatori- 
ally trapped hot ions set up a potential barrier for the 
cold ion beams and facilitate formation of electrostatic 
shocks by reflecting them from the equatorial region. Sim- 
ulations with and without the hot plasma show different 
flow properties; the formation of electrostatic shocks oc- 
cur only in the former case. The simulation with the hot 
plasma also reveals that the magnetic trapping in con- 
junction with the evolution of the electrostatic potential 
barrier produces ion velocity distribution functions con- 
sisting of a cold core and a hot ring in the perpendicular 
velocity. Such a distribution function provides a source 
of free energy for equatorial waves. The corresponding 
electron population is warm and field aligned. 

Introduction 

So far most theoretical studies on plasmaspheric refill- 
ing have been primarily concerned with the outflow of cold 
ionospheric plasma and its trapping in the flux tube. In 
such theoretical studies an important observational fact, 
which has been ignored, is that the flux tubes undergo- 
ing refilling contain a hot plasma population trapped in 
the equatorial region. Such plasmas originate from the 
ring current or the plasma sheet and are characterized by 
Tf[ > 7^, where T*[ and are the perpendicular and 
parallel temperatures of the hot ions, respectively. Spe- 
cific observations of such hot ions having relatively large 
pitch angle anisotropies (A, = T/[/T x ^ > 2) come from 
GEOS-1 and -2, which observed the hot ions in the en- 
ergy range > 10 keV in the noon sector. Such hot ions 
are known to excite electromagnetic ion cyclotron waves [ 
Roux et al., 1982]. 

Another set of observations, where the role of hot 
plasma has been invoked, deals with thermal ions trans- 
versely heated to energies up to a few hundred eV and 
trapped in the equatorial region [Olsen et a/., 1987]. Here 
again hot ions are suggested to be the source of free en- 
ergy for exciting the broadband waves observed from DE- 
1. Theories suggest that the broadband waves are driven 
by a combined effect of temperature anisotropy and ring 
type of distributions of energetic protons [Perraut et ai f 
1982 ]. The effect of heating of the thermal ions on re- 
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filling has been studied [Singh and Chan , 1992], but the 
direct effect of the hot ions on the refilling has not been 
studied so far. 

The purpose of this letter is to show that the hot 
anisotropic ions, which are commonly present in the equa- 
torial region of flux tubes undergoing refilling and drive 
the equatorial processes discussed above [Olsen 
et al , 1987; Roux et al . 1982], facilitate the process of 
electrostatic shock formation since they provide an effec- 
tive potential barrier for the upflowing ion beams of cold 
ionospheric plasma. The shock formation occurs even for 
highly supersonic ion beams which are expected to occur 
during early stage refilling [Banks et al, 1971; Singh et al , 
1986; Rasmussen and Schunk, 1988; Singh , 1990; Wilson 
et al , 1992]. If the hot plasma consists of isotropic elec- 
trons and anisotropic biMaxwellian ions with anisotropy 
^4- = T?JT§ > 1, the potential difference between the 
equator, where the minimum magnetic field strength is 
B m , and an ionospheric point where the magnetic field 
is B, is given by $ = M’jf/efl + T^jT,]- x ln(r), where 
T = A;(l - B m /B) + B m jB [Whipple, 1977], 

Since Bm/B << 1, the equatorial potential with respect 
to the ionosphere for T e << T* is $ ~ (kTJe) ln(A,). For 
Ai — 2 and kT t /e - 10 V, the typical values for these pa- 
rameters used in wave analysis [Roux et al } 1982], $ ^ 
7 V. However, this potential difference is true when no 
ionospheric cold plasma is present; in the presence of a 
cold plasma, it is not certain how large the potential dif- 
ference is and how it is distributed along the field line. 
In the following discussion, we present results from nu- 
merical simulations elucidating the interactions between 
cold and hot plasmas when the former plasma flows into 
the latter one trapped in a magnetic mirror. The simula- 
tion presented here is small scale, in contrast to the large- 
scale problem in space. Therefore, the results presented 
here only serve the purpose of elucidating the processes 
involved in hot-cold plasma interactions. 

Numerical Model 

We perform a one-dimensional particle-in-cell simula- 
tion of plasma flow along an artificial flux tube (Figure 
la). The magnetic field B(X) = B 0 (l — a exp[— (AT — 
d/2) 7 /cr 2 }) where B 0 is a constant field outside 
the minimum-field region, d is the size of the simulation 
system and the choice of ct and <r determines the desired 
field distribution. The hot plasma is created by injecting a 
large number of electrons and ions in the minimum-B field 
region. Such plasma particles are chosen from Maxwellian 
or bi- Maxwellian velocity distributions depending upon 
the requirement of a simulation run. The cold plasma 
flows into the flux tube from the two plasma reservoirs at 
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the ends of the simulation system at X = 0 and X = d 
(Figure la). The simulation technique is described in 
Singh and Chan [1992]. We solve for the motions 
of charged particles in self-consistent electric fields, de- 
termined by solving Poisson’s equation. As the particles 
move in the flux tube, their magnetic moments are as- 
sumed to be conserved. 

In the cold-plasma reservoirs electron and ion temper- 
atures are To. The reservoirs supply a continuous flux of 
charged particles into the flux tube through the process 
of plasma expansion. In the simulations reported here, 
we have used rrii/m, = 400, which adequately separates 
the electron and ion time scales and, at the same time, 
allows computationally feasible runs. The hot plasma is 
injected into the central region of the simulation system 
at time t=0 (Figure la). The properties of the hot plasma 
are described later. Numerical parameters of the simula- 
tions are as follows: system size d = 5 xlO 3 A<*, B field 
parameters a = 0.9, <r = 750 A dy cell size Ax = 20 A d , and 
time step At = O.lu^J,, where Aj and are the Debye 
length and electron- plasma frequency in the cold plasma 
reservoirs. In the following discussion we have used nor- 
malized quantities defined as follows: time t = fuy 0 , dis- 
tance X = X/A d, velocity V = V jV iz and electric poten- 
tial $ = e$/fcT 0 , where V te = ( kTo/m e ) 1 ^ and = 
(m r /m i ) ,/, u> peo . 


Numerical Results 

First we present results from a simulation in which 
“equatorial” hot plasma is not included. This simula- 
tion serves as a reference against which the hot plasma 
effects on the cold plasma flow are compared. Figures lb 
to If show the temporal evolution of the flow of cold ions 
into the flux tube from the reservoirs shown in Figure la. 
These panels show the phase-space plots in the X — Vj| 
plane. Each dot in the figure represents an ion, giving 
its parallel velocity and position. The simulation begins 
at t = 0, when plasmas from the reservoirs begin to flow 
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Fig. 1. (a) Artificial flux tube in which plasma expands 
Rom the cold plasma reservoirs at X=0 and X=d. The 
hot plasma is injected in the central region, (b) to (f) 
ion phase-space plots in X — Vj| plane at selected times as 
shown. The plots show evolution of the cold plasma flow 
in absence of a hot plasma. 


in the flux tube. The corresponding flow of electrons is 
not shown here. Figure lb shows expanding ion beams 
into the simulation region at t = 200. The plot for i = 
1,000 (Figure lc) shows that the ion beams have expanded 
from one end to the other end of the simulation region, 
setting up counterstreaming. Equatorial shock formation 
has not occurred. At this time, ion beams are too fast 
to couple through the ion-ion instability; the relative flow 
velocity between the beams is V rti = 0.4V t e - 8 c # , where 
c i s the local ion-acoustic speed at the equator, X = d/2, 
and is about .0bV te . However, at later times the beams 
become progressively slower as the plasma accumulates in 
the flux tube, and when they become sufficiently slow in 
the regions near the ends of the flux tube they excite ion- 
ion instability as seen for t > 1,400 (Figures Id to If); 
the instability creates vortex-like structures in the phase 
space. 

The instability is seen to thermalize the ion beams (Fig- 
ure le). However, in the middle of the simulation region, 
where the magnetic field is minimum, counterstreaming 
of ion beams is seen to persist. 

Simulation with a hot plasma in the region of mini- 
mum magnetic field reveals a quite different behavior of 
the flow. We performed several simulations by varying the 
temperatures of the hot plasma but kept ion anisotropy 
A, >1 and electron anisotropy A t = 1 and T t « T*. 
Such properties of hot plasma have been measured [Roux 
tt al , 1982]. Relative densities of the cold (n c ) and hot 
(ti#) ions in their respective source regions were also var- 
ied. The results from such a parametric study will be re- 
ported elsewhere. We present results here for 7^ = T^/2 
= 900T o , T e = 10T o , and n c ^ 100 n H - 

The hot plasma is injected in the central region 2300 < 
x/A d < 2700 at t = 0, when cold plasma flows begin from 
the reservoirs. Figure 2 shows the evolution of the cold 
plasma flow along with the potential and density profiles 
in the flux tube. The panels a to d show the ion phase 
space plot in X — Vjj plane. The phase-space plot of ions in 
X — Vj_ plane is shown in panels e to h. The evolution of 
the potential profiles in the flux tube is shown in panels i to 
I, and the corresponding plasma density is shown in panels 
m to p. The evolution of the electron phase-space (X - Vjj) 
is shown in panels q to t. As shown on the top of the figure, 
the columns from left to right show the properties of the 
flow at f = 100, 500, 700, and 1000, respectively. At an 
early time (t = 100), the flow can be characterized by the 
following features: (1) expanding cold ion beams (panel 
a), (2) equatorially trapped ions originating from the hot 
plasma (panel e), (3) a potential maximum at the mid- 
point of the simulation region (panel i), and (4) a density 
maximum coinciding with the potential maximum (panel 
m). The density maximum at this time entirely represents 
the trapped hot ions, and there is no contribution from the 
cold plasma flow. 

At t = 500, the cold plasma flow comes into contact 
with the hot plasma; the cold ion beams have somewhat 
penetrated into the region of the hot plasma (panels b 
and f). Already the cold beams show the sign of retar- 
dation by the potential barrier. The slowing down of the 
ion beams enhances the potential barrier further; panel j 
shows that at t = 500 the equatorial potential maximum 
has grown to 15fcTo/e, in contrast to 3kTo/c at t = 100. 
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Fig. 2. Evolution of the cold plasma flow in the presence of equatorially trapped hot ions. Each column shows the 
state of the plasma at times shown at the top. (a) to (d) -Y — Vjj phase-space plots for ions, (e) to (h) X V± phase 
space plots for ion, (i) to (1) Potential distributions, (m) to (p) Ion density distributions, and (q) to (t) Electron 

phase- space plots in .Y — Vjj plane. 

propagation of shocks can be seen by comparing the pan- 
els in columns for t = 700 and 1000. The shocks propa- 
gate with a speed V $h ~ 0.09V,. ~ C to , where C to is the 
ion-acoustic speed in the cold plasma reservoirs. As the 
shocks propagate away from the hot plasma, the thermal- 
ized cold plasma behind the shocks punches through the 
region of hot plasma where the potential barrier existed 
earlier. This penetration of the cold plasma is clearly 
seen from the X - V L plots shown in Figure 2e to 2h. 
As the cold ions penetrate into the minimuni-B field re- 
gion, where hot ions are trapped, they cool adiabaticily 
as clearly seen from these plots. The penetration of cold 
plasma is accompanied with a reduction in the potential 
barrier, caused by the cold electrons being accelerated into 
the high potential region (panels r and s); these electrons 
tend to neutralize the effect of the hot ions stably trapped 
in the minimum-B field region. 

A noteworthy feature of the plasma after the spatial 
mixing of hot and cold plasmas is the nature of the re- 
sulting ion distribution functions in the equatorial region, 
as shown in Figures 3a and 3b. The parallel velocity dis- 
tribution function (Figure 3a) shows a core of cold ions 
superimposed on a warm ion population. The perpendic- 
ular velocity distribution function (Figure 3b) also shows 
the cold ions, but the hot ions appear as a beam. Since 
the ions in the beam are nearly uniformly distributed in 
their phase, the beam is actually a ring in the perpendic- 
ular velocity space. The ring distribution function is the 
result of the trapping of the hot ions in the magnetic mir- 
ror in combination with the evolution of the electrostatic 
potential distribution. Some of the ions from the hot pop- 
ulation having relatively small perpendicular velocities are 
lost from the mirror due to the parallel electric fields. In 
the absence of the electrostatic potential, one expects a 


Note the change in the vertical scale between panels i and 
j. The filling of the flux-tube with the cold plasma changes 
the density profile by creating a density minimum at the 
equator (panel n). 

The plots at i = 700 show that the contact between 
the hot and cold plasma has evolved into a pair of elec- 
trostatic shocks, one on each side of the hot ions trapped 
in the minimum-B field region. The shocks are indicated 
by arrows where flow velocities of the cold ion beams sud- 
denly decrease (panel c). Near the shocks, the potential 
profile shows sharp jumps (panel k) and the density profile 
shows spikes (panel o). 

The shocks propagate away from the “equatorially” 
trapped hot plasma towards the end of the flux tube, 
thermalizing the cold plasma just behind the shocks. The 
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Fig. 3. Ion velocity distribution functions in the equato- 
rial region (a) Parallel velocity distribution, (b) Perpen- 
dicular velocity distribution. 
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loss-cone distribution. The loss of ions with relatively low 
yields the ring distribution. The cold core originates 
from the cold ions penetrating in the minimum-B field re- 
gion. The importance of these findings lies in the fact that 
ring types of distribution functions along with a core of 
cold ions can be a source of free energy for exciting equa- 
torial waves [Lee and Birdsall , 1979], which can heat the 
latter ions. We point out that the low energy end of the 
ring distribution found here occurs at an energy of a few 
tens of eV, unlike the relatively energetic ring distribution 
with energy > 5 keV observed from GEOS [Perraut et al., 
1982]. 

Conclusion and Discussion 

The main conclusions of the study in this paper are the 
following: (1) The presence of the hot anisotropic ions 

trapped in the equatorial region can drastically affect the 
flow of the cold plasma. (2) The potential barrier associ- 
ated with such a hot ion population facilitates formation 
of electrostatic shocks even when the cold ion beams com- 
ing from the ionosphere are highly supersonic. (3) The 
eventual mixing of the hot and cold plasmas produces a 
potentially unstable velocity distribution function for the 
ions. They develop a ring distribution in the perpendic- 
ular velocity. The ring distribution with df(V±)/dV± be- 
gins at a relatively low energy corresponding to the energy 
of an ion falling down the potential barrier of the order of 
10 V. The ring distribution discussed here is different from 
the observed high energy rings (> 5 keV) produced by hot 
plasma injection and its subsequent convection [ Perraut et 
a/., 1982]. 

Recent satellite observations have shown that equato- 
rially trapped warm ions are generated by equatorial per- 
pendicular heating of the thermal ions; the heating is 
caused by waves generated by an energetic hot ion popula- 
tion, having the temperature anisotropy T * |[ > and/or 
ring type of velocity distribution [Olsen et al, 1987; Per - 
raut et a/., 1982 ]. Olsen et al. [1987] also report the 
reflection of cold ion beams. However, it is believed that 
the reflection is caused by the potential barrier set up 
by the thermal ions which have undergone the transverse 
heating. From our simulations, it appears that the hot 
ion population, which drives the waves by virtue of its 
anisotropy or ring, may reflect the fast ion beams at an 
earlier stage before the waves grow and consequent ion 
heating produces a warm anisotropic ion population. 

Futhermore, before thermal ions are heated they must 
enter the hot plasma region. This entry is retarded by the 
potential barrier set up by the hot ions. Even if the hot 
plasma injection occurs at a stage when field-aligned flows 
have set up at the equator, the potential barrier is likely 
to expel the cold plasma from the equatorial region. The 
presence of thermal ions is essential in the wave generation 
processes. Therefore, it is suggested that the reflection of 
ion beams observed from DE-1 may not be necessarily 
caused by the equatorial perpendicular heating of ther- 
mal ions; it appears that the reflection of ion beams and 
heating of thermal ions are all driven by the hot ions, but 


heating should phenomenologically follow the reflection of 

upflowing ion beams. 
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ABSTRACT 


Hydrodynamic and semi-kinetic treatments of plasma flow along closed geomagnetic 
fields lines are compared. The hydrodynamic treatment is based on a simplified 16- 
moment set of transport equations as the equations for the heat flows are not solved; 
the heat flows are treated heuristically. The semi-kinetic treatment is based on a 
particle code. The comparison deals with the distributions of the plasma density, 
flow velocity, and parallel and perpendicular temperatures as obtained from the two 
treatments during the various stages of the flow. In the kinetic treatment, the appro- 
priate boundary condition is the prescription of the velocity distribution functions 
for the particles entering the flux tubes at the ionospheric boundaries; those parti- 
cles leaving the system are determined by the processes occurring in the flux tube. 
The prescribed distributions are half-Maxwellian with temperature To and density 
no* In the hydrodynamic model, the prescribed boundary conditions are on den- 
sity (n 0 ), flow velocity (V 0 ) and temperature (T 0 ). We found that results from the 
hydrodynamic treatment critically depend on Vo; for early stages of the flow this 
treatment yields results in good agreement with those from the kinetic treatment, 
when Vo = yj kTo/2irm, which is the average velocity of particles moving in a given 
direction for a Maxwellian distribution. During this early stage, the flows developing 
form the conjugate ionospheres show some distinct transitions. For the first hour or 
so, the flows are highly supersonic and penetrate deep into the opposite hemispheres, 
and both hydrodynamics and kinetic treatments yield almost similar features. It 
is found that during this period heatflow effects are negligibly small. When a flow 
penetrates deep into the opposite hemisphere, the kinetic treatment predicts reflec- 
tion and setting up of counterstreaming. In contrast, the hydrodynamic treatment 
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yields a shock in the flow. The reasons for this difference in the two treatments is 
discussed, showing that in view of the relatively warm ions, the coupling of ion beams 
and the consequent shock formation in the offequatorial region are not likely due to 
the enhancements in the beam temperatures. The counterstreaming in the kinetic 
treatment and the shock in the hydrodynamic treatment first advance upward to 
the equator and then downward to the ionospheric boundary from where the flow 
originated. The transit time for this advancement is found to be about 1 hour for 
the respective models. After 2 hours or so, both models predict that the flows from 
the ionospheric boundaries are generally subsonic with respect to the local ion-sound 
speed. At late stages of the flow, when a substantial fraction of ions entering the 
flux tube begin to return back in the kinetic treatment, the hydrodynamic treatment 
with the boundary condition Vo = ^kToj^hrm yields an over-refilling, and the choice 
of Vq becomes uncertain. 
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1. Introduction 


In connection with the problem of plasmaspheric refilling, in recent years several 
models for plasma flow along closed magnetic field lines have been developed [Khaz- 
anov et al., 1984; Singh, et al., 1986; Singh , 1988; Rasmussen and Schunk , 1988; 
Singh, 1991; Guiter an d Gombosi, 1990; Wilson , et al , 1992]. These models dif- 
fer in complexity in terms of describing the plasma and in including the ionosphere 
as a source of plasma for the refilling. In terms of describing the plasma the most 
contrasting feature of the existing models deals with the hydrodynamic and kinetic 
treatments for the flows, based on plasma fluid equations and a particle-in-cell code, 
respectively. For the purpose of including the ionosphere as a source of plasma for 
the refilling, in most studies the topside ionosphere is replaced by a set of bound- 
ary conditions on the plasma flow, except in Guiter and Gombosi [1990] who have 
included the generation and loss of plasma through chemical reactions in a hydro- 
dynamic model. In this paper, we are mainly concerned with the treatments of the 
plasma with simple sets of boundary conditions on the plasma flow; we compare 
the properties of the plasma flow along closed field lines as given by hydrodynamic 
[Singh, 1992] and semi-kinetic [Wilson et al., 1992] treatments. 

It is generally believed that the kinetic treatment of a plasma is superior to a 
hydrodynamic (fluid) one. The success of a hydrodynamic treatment depends on the 
problem being solved and on the ingenuity of the user in chosing the hierarchy of 
moment equations on which fluid equations are based. In recent years, researchers 
in space physics have used fluid descriptions based on 13-moment [ Schunk , 1977; 
Mitchell and palmadesso, 1983], and 16-moment [Barakat and Schunk , 1982; Gan- 
guli and Palma c itsso , 1987; Gombosi and Rasmussen, 1991; Rorosmezey et al., 1992, 
1993] set of transport equations. In developing the moment equations, the ingenuity 
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lies in a series expansion of the plasma distribution function using a biMaxwellian dis- 
tribution function as a base. Therefore, hydrodynamic treatment based on moment 
equations are good as long as the distribution function is close to a biMaxwellian. 

9 

When the distribution function severely departs from a biMaxwellian and involves 
multistreaming of plasma particles, the moment equations are seriously handicapped, 
despite the sophistication incorporated through the use of higher order moment equa- 
tions. 

As mentioned above, recent models for plasmaspheric refilling are based on both 
a kinetic treatment using PIC code and hydrodynamic treatment with varying de- 
gree of sophistication in chosing the hierarchy of the moment equations. In some 
early models only continuity and momentum equations were solved for the ions and 
the electrons were assumed to remain isothermal [Singh et al. , 1986; Rasmussen and 
Schunk , 1988; Singh , 1988]. Studies which included temperature equations assumed 
that either the heatflow is given by the collision-dominated thermal conductively 
[Khazanov et al 1984; Guiter and Gombosi , 1990] or ignored the heatflow com- 
pletely [Singh, 1992; Singh and Chan , 1992]. Neither of these assumptions correctly 
describe the heat transport in the refilling problem [ Singh and Horwitz , 1992]. In 
the collisionless limit of plasma flow during refilling, the usual description of heat- 
flow in terms of Spitzer thermal conductivity breaks down and such a treatment 
overestimates the heatflow. When the heatflow is handled by the 16-moment set of 
equations, the problem remains in the sense that the heatflow equations are valid for 
only relatively small heatflows. When the heatflow becomes large, the validity of such 
equations ceases and numerical instabilities result in computational work. Since, a 
priory it is not known when a large heatflow develops in a model, ad-hoc damping 
mechanisms are included to damp the heatflow, whether it is physically warranted 
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or not [Palmadesso et al. 9 1988; Rasmussen, Private Communications]. The above 
problem with the 16-moment transport equations is generic; it arises irrespective of 
the sophistication in numerical techniques employed in solving them [ Korosmezey et 
a/., 1992, 1993). 

Despite the above difficulties with the hydrodynamic treatment, it has been used 
for practical reasons because it provides simplicity and considerable economy in com- 
putational work and depending on the plasma conditions it can work successfully. 
Therefore, it is advisable to keep in mind the assumptions made in using this treat- 
ment and, if possible it is even better to check the validity of this treatment by 
comparing its prediction against that from a kinetic treatment. Such a comparison 
may reveal when and how a fluid model succeeds. 

The purpose of this paper is to carry out a comparison between models of the 
plasma flows along closed field lines based on kinetic and hydrodynamic treatments. 
The former treatment uses a particle-in-cell (PIC) code for ions [ Wilson et al , 1992]. 
The latter one uses transport equations for the flow of mass, momentum, and parallel 
and perpendicular temperatures of ions [Singh, 1992], but the heatflow is treated 
heuristically [Metzler, 1979]. In both the treatments electrons are assumed to obey 
the Boltzmann law. In the present paper, the ionospheric outflows is included by 
imposing a set of boundary conditions on the flow of ions at an altitude of 2000 km. 
The choice of this altitude is primarily due to the existing models [Wilson et a/., 
1992; Singh . 1992] in which ionospheric loss and generations processes for the plasma 
are not yet included. 

The closed field lines provide the possibility of a variety of flow conditions ranging 
from highly supersonic to subsonic flows as an empty flux tube refills. Futhermore, 
the flows along closed fields fines develop counterstreaming due to interhemispheric 
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plasma flows. Since hydrodynamic treatments are most suspect under counterstream- 
ing situations [ Manheimer et a/., 1976], the comparison carried out here provides a 
useful guide for assessing the validity and usefulness of a hydrodynamic treatment. 

We have found that for the conditions of highly supersonic flows, the two-stream 
hydrodynamic treatment yields flow properties in good agreement with that from 
the semi-kinetic treatment. Demars and Schunk[1991] reported a similar agreement 
based on 16-moment set of equations including heat flows. We find that the bulk 
parameters such as the density, flow velocity and temperatures are in good agreement 
even for the simpler hydrodynamic model when heatflow is included heuristically; the 
reason being simply that when the flow is highly supersonic, the dominant transport 
of heat is through the bulk flow velocity and the transport due to the thermal effects 
is negligibly small. 

When reflection of flows causes counterstreanung, the hydrodynamic treatment 
gives rise to shock formation, which is not seen from the kinetic treatment. How- 
ever, when the counterstreaming flows become subsonic, the hydrodynamic and semi- 
kinetic treatment again produce flow properties in reasonable agreement. 

In order to study plasma flow in space, the plasma treatment must be supple- 
mented by a set. of boundary conditions on the flow equation. For the closed field lines, 
the boundary conditions are determined by the top-side ionosphere. The boundary 
conditions along with the demand for plasma at high altitudes produce the flow. The 
ionospheric boundary conditions involves generation and loss of ionospheric plasma 
particle species. Since here our primary goal is in identifying the kinetic and fluid-like 
behaviors of plasma flow and not the supply of plasma from the ionosphere and the 
refilling rate, we have simulated the outflow of ionospheric plasma by imposing a set 
of boundary conditions at an altitude of 2000 km in both the hemispheres. In the 
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semi-kinetic model, the imposed boundary condition is on the velocity distribution 
function of the ions entering the flux tube. It is assumed to be half- Maxwellian. 
The returning particles are self-consistently determined. In the hydrodynamic model 
the boundary conditions are the moments of such a distribution. Since the kinetic 
effect dealing with the returning particles are lost in the hydrodynamic model, the 
hydrodynamic model does not agree with the kinetic model when returning ion flux 
becomes sufficiently large. Can this disagreement be resolved by a more sophisticated 
treatment of the heatflow by using a complete set of 16-moment equations? In order 
to answer this question, further comparative studies are suggested. 

The rest of the paper is organized as follows. The theoretical models are described 
in section 2. The comparison between the results from the two models is carried out 
in section 3. The main conclusions of the paper and their discussion are given in 
section 4. 


( 


2. Theoretical models: 


The semi-kinetic model, which is based on a particle-in-cell code, has been pre- 
viously described for both open [ Wilson , et al 1990] and closed [ Wilson, et ai, 
1992] flux tubes. Coulomb collisions are included in the model, the collisions are 
implemented by pairing simulation ions' according to an algorithm which conserves 
energy and momentum [Takizuka and Abe, 1977]. The algorithm yields good ap- 
proximation for the collisions when the collisional relaxation time is shorter than the 
time step in advancing the ion motion. In the hydrodynamic model, we solve the 
plasma transport equations based on 16-moment approximation [e.g., see Ganguli 
and Palmadesso , 1987, and Barakat and Schunk , 1982] 
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where t is time; r is the geocentric distance along the flux tube; s is the distance 
along the tube from the northern ionospheric boundary at A = + \o (see Figure 1), 
n, V, Tj| and Xj_ are the number density, flow velocity and parallel and perpendicular 
temperatures of ions in the plasma flow, respectively; 
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(jll and q± are the heat fluxes along the magnetic field line associated with Tj| and 
Tx, respectively; E is the parallel electric field; gf|| is the component of the gravita- 
tional force parallel to the magnetic field, and rrt and e are the ion mass and charge, 
respectively. The collision terms denoted by [*] c are calculated using Burger s Formu- 
lae [Burger, 1969], which are modified to include flow velocity corrections [Mitchell 
and Palmadesso et al. , 1983; Ganguli and P almadesso , 1987] and the correction for 
temperature anisotropy [Ichimaru et a/., 1973; Singh , 1991]. 

We do not solve the heat flow equations, which have proven to be quite trou- 
blesome to solve numerically [ Palmadesso , et ai } 1988]. The difficulty arise for a 
relatively large heatflow, for which the moment equations themselves become invalid. 
Since it is unpredictable in a model when the heatflow may be large, ad-hoc pro- 
cedures are employed to attenuate the heatflow for the numerical stability of the 
models. This has been found to be true irrespective of the numerical techniques 
for solving the equations [Palmadesso et a/., 1988; Roromezey et al. y 1992, 1993; 
Rasmussen, private communication]. 

In this paper instead we have included the effects of heat flow heuristically by 
closely following the treatments in solar wind studies [e.g. Metzler , et al 19 1 9] . In 
a collisionless plasma the usual picture of heat flow, given by q a — -A a Vr a with 
K a as the thermal conductivity, may not be valid because L, the mean free path, 
is > L t = (T- l dT/ds)-\ the scale length in the temperature variation. In such a 
collisionless situation, the heat flux can be calculated on physical ground as follows. 
The heat fluxes and across a surface in a single direction in a plasma described 
by a biMaxwellian distribution function with parallel and perpendicular temperatures 
7j| and 7T, respectively, say along the magnetic field vector, are given by 

qjf = nJfcT„(fc7]|/27rm) 1/2 (5) 
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q u t = nkT i_(kT\\j 2nm) 


( 6 ) 


1/2 

In a uniform plasma for which the distribution function is independent of the parallel 
coordinates, the heatfiow at any point is zero because heat flux in a given direction is 
cancelled by the heat flux in the opposite direction. In the presence of a spatial inho- 
mogeniety, the cancellation is not complete, and the heat fluxes q\\ and appearing 
in equations ( 3 ) and ( 4 ) can be heuristically written as [ Metzler , et al., 1979], 

q a ~ er]nkT C 'Vt\\, (7) 

where the subscript a stands for J_ or |], e = — 1 if dT a /ds > 0 and t = 1 if 
dT a /ds < 0 . Thus, in the heat flow model adopted here only the sign of the heat flux 
depends on the temperature gradient and not its magnitude. The factor tj determines 
the reduction in the heat flow below the unidirectional fluxes in (5) and ( 6 ). Later on 
in this paper we show that 77 in the range say 0.1 - 0.3 yields results in a reasonable 
agreement with the kinetic model, in which heat fluxes appear self-consistentiy. A 
similar model for heatfiow was used by Singh[1992] for plasma flow along open field 
lines. In both the hydrodynamic and kinetic models adopted here, electric field E is 
calculated by assuming that the electrons obey the Boltzmann law and the condition 
of quasi-neutrality prevails. 

The plasma flow along a closed field line is studied by solving an initial-boundary 
value problem. In the hydrodynamic model, the plasma flows originating from 
the conjugate ionospheres are treated as separate fluids; this treatment is termed 
as a two-stream model [Singh, 1988; Rasmussen and Schunk , 1988; Singh , 1990]. 
In both the models, it is assumed that at the initial time (t = 0) the flux tube 
is highly depleted. The depletion is given by n t = raisin A/ sin A 0 ) w i t h the 
minimum density limited to 10~ 4 n 0 , where ng is the density at the ionospheric 
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base A = ±A 0 (Figure 1). Initial flow velocity V(X ,f = 0) = 0 and temperatures 
T±{\,t = 0) = = 0) = T 0 = 0.3 eV. In the hydrodynamic model, the 

boundary condition for the fluids originating from the northern hemisphere are. 
n n (X = A 0 , 0 = n 0 ,V n (X = X 0 ,t) = T„ n (A*, f) = T^A,,*) = T 0 \ at the bound- 

ary A = —X 0 floating boundary condition are applied. A set of similar boundary 
conditions are used for the fluid originating from the southern hemisphere, but with 
the roles of A = ±A 0 interchanged. In the kinetic model, the boundary conditions on 
ion distribution function /(A, V) are that f(X = X 0 ,V>0) and f{X = —X o ,V<0) 
are half-maxwellians, with a temperature To. These boundary conditions prescribe 
only the ions entering the flux tubes. The ions leaving the flux tubes are deter- 
mined by the processes occurring inside it. A half-Maxwellian, and not a displaced 
Maxwellian, is chosen because of the following reasons. There is no clear observa- 
tional evidence of Supersonic flows along closed field lines at an altitude of 2000 km. 
Futhermore, our calculations show that in about 2 hours the flow in the flux tube 
becomes subsonic nearly all along the flux tube; only for an initial stage of about 2 
hours, supersonic flows are seen. In view of such uncertainties on the the flow velocity 
at 2000 km. a half-Maxwellian serves the purpose of the comparative study. 
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3. Numerical Results 


We compare here the properties of the flow in a flux tube with L = 4 as seen 
from the semi- kinetic and hydrodynamic models. Since the boundary value of the 
flow velocity (Vo) and the heat flow factor 77 in the hydrodynamic model are free 
parameters, comparison is performed by varying them over physically reasonable 
ranges. The comparison also deals with the accumulation of plasma in the flux tube 
and the equatorial plasma density. 

3.1 Initial Supersonic Flow 

We recall that the hydrodynamic model is based on two-stream flow in which flows 
originating form the two hemispheres are treated as separate fluids, and the temporal 
evolution of the two streams is separately studied. Likewise, even in the semi-kinetic 
model, the separate identity of the ions originating from the two hemispheres is 
maintained. Figure 2 shows the evolution of the flow originating from the northern 
hemisphere as seen from the semi-kinetic model. This figure gives the phase-space 
density plots in A — V]| plane, where A is the geomagnetic latitude and Vj| is the 
flow velocity along the magnetic field line. The positive and negative values of A 
correspond to the northern and southern hemispheres, respectively. The darkest 
region in the grey-scale plots represents the highest density as indicated by the scale 
on the right-hand side. At time t = . 003 hour , the plasma in the tube is essentially the 
initial plasma with a density profile given by n = raisin A/ sin \o} 6 • At later times this 
plasma expands into the flux tube and it is seen to cross the equator at t = 0.25 hour. 
Along with the expansion, new plasma enters the flux tube at the boundary A = -fAfl. 
It is seen that by the time t = 0.75 hour , the flow has penetrated all the way to the 
opposite boundary at A = — \ 0 . It is found that the plasma reaching this boundary 
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is not totally lost, but it is partially reflected back, setting up a counterstreaming 
flow as seen from the plots for t > 1 hour. The reflected flow is seen to reach the 
boundary at A = \ 0 by the time t = 2 hours. The plasma flow originating from the 
southern hemisphere shows a similar behavior as shown in figure 2, with the role of 
boundaries at A = ± Ao interchanged. It is worth pointing out that after reflections, 
ions merge with the ion stream moving in the opposite direction and they do not 
appear as a separate ion beam. 

The possible consequences of the counterstreaming flow will be discussed later 
on. We now compare the above features of the flow seen from the kinetic model with 
those seen from the hydrodynamic model. Figures 3a to 3d show the comparison for 
t = 0.5 hour; these figures show the distributions of (a) density, (b) flow velocity, (c) 
parallel temperature, and (d) perpendicular temperature. In each panel the curve 
from the kinetic model is labeled, and the curves from the hydrodynamic model 
for three values of the heatflow reduction factor tj = 0 , T] = 0 .05 and r/ = 0.3 are 
indicated by the legend. It is seen from these figures that for most of the flux tube all 
four curves are quite close together, irrespective of the heatflow factor 77 . However, 
near the opposite boundary (A = — A 0 ), the curves from the hydrodynamic model 
tend to diverge from the kinetic model, depending on the value of 77 . This difference 
between the two models is attributable to the fact when the flow begins to slow down 
due to a relative increase in the plasma density, the hydrodynamic model predicts 
an increase in the parallel temperature as clearly seen for A < —36° in Figure 3c. 
The increase in Tji enhances the pressure and further slows the flow and enhances 
the density. In the semi-kinetic model, the temperature enhancement does not occur. 
Instead, a counterstreaming develops. This contrast between the two models becomes 
much clearer at later times, for example, at t = 1 hour for which the comparison is 
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shown in Figures 4a to 4d. 

Figure 2 shows that at t = 1 hour the reflected ions set up counterstreaming 
throughout the southern hemisphere (A < 0). Since the hydrodynamic model cannot 
handle the counterstreaming, the reflection process creates a shock, which is clearly 
seen in the density, velocity and temperature plots in Figures 4a to 4c, respectively, 
across the shock indicated by the arrows, density suddenly increases, flow velocity de- 
creases and the parallel temperature also increases. We note that the hydrodynamic 
curves for different values of rj begin to show some difference among themselves, 
with the curve for 77 = 0.3 being closest to that from the kinetic model. It is worth 
pointing out that the flow velocity in the kinetic model is the average over the coun- 
terstreaming ions. The average velocity is lower than that from the hydrodynamic 
model over the region of the counterstreaming, but where the counterstreaming has 
not yet occurred (A > 30°) the flow velocity from the two models are generally in 
good agreement. 

The shock first propagates upward to the equator and then downward and reaches 
the ionospheric boundary at A = \o at t ~ 2 hours. The propagation of the shock 
in the density profile is shown by the arrows in Figure 5. The transit time of about 
2 hours for the shock is in agreement with the development of the counterstream- 
ing starting in the southern hemisphere and spreading to the northern ionospheric 
boundary by t = 2 hours, (see Figure 2). We find that the heatflow plays only a 
minor role in the motion of the shock; the shock speed is slightly enhanced with 
increased heatflow 77; for 77 = 0.3 and 0.05 the shock is already absorbed near the 
boundary A 2^ A^ , while for 77 = 0, the shock can be still seen in the flux tube at t — 
2 hours . 

After the shock reaches the boundary at A = A 0 , the flow in the flux tube becomes 
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generally subsonic with respect to the ion- acoustic speed, which is about 10 km/s 
with electron and ion temperatures T 0 = 0.3ev at A = ±A 0 - We will discuss the 
subsonic stage after we examine the reason why a shock did not form during the 
early stage ( ~ 1 hour) of the counterstreaming (Fig. 2) in the semi-kinetic model. 

3.2 Electrostatic shock 

We have just seen that a shock automatically forms in the hydrodynamic model 
as soon as the flow begins to reflect near the opposite boundary. On the other hand, 
the kinetic model does not show the shock formation. Instead, a counterstreaming 
flow develops. We examine this issue in terms of the conditions for shock formation 
and ion velocity distribution function. 

According to the original suggestion of Banks et a/.[19/l], a shock should form 
when supersonic flow from the conjugate hemispheres collide at the equator. The 
flows collide as early as t = 0.25 hour ; the flows from the northern hemisphere can 
be seen from Figure 2 and the corresponding flow from the southern hemisphere is 
the mirror image of this flow with respect to the equator. When the flow begins to 
overlap, the shock should form through the ion-ion instability. The conditions for 
such an instability in a colliding situation are given by 

1 -3v tt < Vb<MC, and T t >3T\\ (8) 

where V ft is the ion thermal velocity, V b is the ion beam velocity, C s is the ion- 
acoustic speed and M is the mach number, which could be as large as 4 [Forslund 
and Shonk , 1970; Montgomery and Joyce , 1969]. However, it must be mentioned here 
that high critical Mach number M Cr J is determined by the non-linear evolution of 
the electron dynamics including trapping and heating of electrons [e.g., see Singh , 
1988]. For isothermal electrons, as it is assumed in the semi-kinetic model, M = 1 .6 
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[Tidman and Krall , 1971]. It is worth pointing out that for large beam velocities, 
oblique ion waves propagating at an angle from the magnetic field are likely to be 
excited. However, the role of such highly oblique waves, which are likely to occur for 
highly supersonic beams, in momentum exchange between interpenetrating beams 
and shock formation is not well understood. 

We examine here the likelihood of the instability occurring from the flow param- 
eters given by the semikinetic model. First we do this exercise for t = 30 minutes 
when the flow has crossed the equator. Figure 6 shows the average flow velocity V*,, 
the temperature ratio Tj|/T e and the ion-acoustic speed C, as function of geomagnetic 
latitude for the flow at t = 30 minutes, shown in Figure 2. Note that the tempera- 
ture ratio is plotted after multiplying it by 10, so that all the plots in Figure 6 can 
utilize the same vertical scale. C, is calculated from C 9 — [k(T e + ZT\\)/m} 1 ^ 2 . The 
critical temperature ratio, Tj|/!T e = 0.33, for the instability is shown by the segment 
of the thick horizontal line in Figure 6. It is seen that the ions have sufficiently 
cooled down to meet instability condition on the ion temperature over an extended 
equatorial region (|A| < 20°). The flow coming from the opposite hemisphere shows 
a similar feature. The ion- acoustic speed in the equatorial region is about 8 km/s. 
It is seen that over the latitudinal region |A| < 20°, the ion beam velocity is about 
Vb ~ 2(7,. In the semikinetic and hydrodynamic models discussed here the electrons 
are assumed to obey the Boltzmann law. Therefore, ion beams with such velocities 
are too fast to excite the ion-ion instability and thereby to form shocks in the model. 
Futhermore, it is important to point out that the processes which lead to shock for- 
mation, including the ion-ion instability, are microprocesses, which are suppressed in 
the large-scale models [Singh and Chan , 1993]. If electrons dynamics were rigorously 
included in the model and the associated microprocesses properly resolved, it is likely 
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that the ion-ion interaction would have occurred forming shocks. 

Figure 7 shows the drift velocity \\ , Cs and T\\/T e aX i = 1 hour for the flow 
originating from the northern hemisphere. It is seen that as the ion beam penetrates 
into the opposite hemisphere (A < 0), it gets progressively warmer and the temper- 
ature condition T\\/T e < 0.3 is not met beyond |A| = 10°. Thus ion instability and 
shock formation are not expected. This indicates that the shock formation in the 
hydrodynamic model (Figures 4 and 5) is an artifact of the model. 

In view of the above discussion in connection with Figures 6 and 7, it emerges 
that on the basis of temperature condition alone, it can be argued that if shocks form, 
they should be during the early stage when the ion beams begin to interpenetrate 
in the equatorial region. During later stages, when the beams penetrate into the 
opposite hemispheres the shock formation is not likely unless some how electrons are 
heated enhancing the temperature ratio T c /Xj|. However, as mentioned earlier the 
shock formation in the equatorial region requires a rigorous treatment of electron 
dynamics. In view of the simplified treatment of electrons in the models described 
here, and relatively coarse spatial and temporal resolutions afforded by them the 
issue of equatorial shock formation can not be settled in this paper. 

i 

3.3 Subsonic flow 

After the initial stage of supersonic flows from the conjugate ionospheres, the 
flows become generally subsonic. This is predicted from both the models . Figure 
8 shows the status of the flow at t = j hours, from both the hydrodynamic and 
semi-kinetic models. As before, there are three curves from the hydrodynamic model 
which are compared against the curve from the semi-kinetic model. Figure 8b shows 
that the flow velocities obtained for different values of rj from the hydrodynamic 
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model agree with the flow velocity given by the kinetic model. The maximum flow 
velocity of about 5 km/s seen near the boundary A = Xo is subsonic with respect to 
the ion-acoustic speed C, = 10 km/s. We note that the average flow velocity peaks 
slightly above the boundary in both hydrodynamic and kinetic models. The peaking 
is a consequence of the boundary conditions at A = Aq and the acceleration of ions 
by the pressure and electric field distributions in the close vicinity of the northern 
boundary of the flux tube. 

Figures 8 a. 8 c, and 8 d show that for the subsonic flow the density and temperature 
structures critically depend on the heatflow factor 77 . For 77 = 0.3, the hydrodynamic 
model yields results in good agreement with those from the kinetic model. When 
77 becomes too small (77 < 0.05), the structures in the density and temperature pro- 
files markedly differ from the kinetic model; the density structure shows an extended 
density cavity in the equatorial region, where parallel temperature is relatively high 
[Singh, 1991]. Futhermore, for low values of 77 there is density enhancement and cor- 
respondingly a low parallel temperature in the southern hemisphere. When heatflow 
factor is sufficiently large (77 > 0.15), such structures in rc(A) and T|j(A) are washed 
away. In a recent paper. Ho et al [ 1993] computed the parameter 77 from the semi- 
kinetic model for flows along open flux tubes and it was found to be > 0 . 1 . 

The comparison between the hydrodynamic and kinetic results at t = 12 hours is 
shown in Figures 9 a to 9 d. The density and temperature structures at this stage are 
qualitatively similar to that at t = 4 hours as shown in Figures 8 a to 8 d. However it 
is seen that at t = 12 hours , the density and temperature profiles even for 77 = 0.05 
have begun to compare well with that for 77 = 0 . 3 , for which the density distribution 
agrees well with that given by the semikinetic treatment. The discrepancy between 
the densities predicted by the kinetic treatment and the hydrodynamic one for 77 = 
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0.3 is bounded by 15 %, for most part of the flux tube, except near the southern 
boundary A — -A 0 . 

The runs for the comparison between the hydrodynamic and kinetic models were 
carried on until t = 4 8 hours. For time t > 12 hours, it was found that the hydrody- 
namic model systematically yields densities higher than that given by the semi-kinetic 
model. Figure 10 shows the comparison between the equatorial densities as obtained 
from the two models . This figure shows the temporal evolution of the equatorial den- 
sities found from the kinetic (solid line) and the hydrodynamic (broken line curves) 
models. For the latter model, the densities are plotted for different values of the 
flow velocity ( Vo ) at the boundaries A = i Ao . We remind ourselves that the results 
from the hydrodynamic model shown in Figures 3, 4, 5, and 8 are for a flow velocity 
V 0 = {kT 0 /2irm) 1/2 = 0J9V t . We notice from Figure 10 that for this boundary 
value of the flow velocity, the kinetic and hydrodynamic curves are remarkably close 
for t < 12 hours. This implies that this boundary value of the flow velocity closely 
corresponds to the input flux determined by a half- Maxwellian distribution function, 
which is imposed as boundary condition in the kinetic model. For t > 12 hours, the 
boundary value of V Q — 0 .39V t yields an over-refilling compared to the kinetic model. 
This simply implies that the net influx of ion into the flux tube at the ionospheric 
boundaries steadily decreases in the kinetic model, primarily due to the ions flowing 
out of the flux tube. On the other hand, in the hydrodynamic model, the influx is 
primarily determined by the imposed flow velocity and it remains constant. This is 
demonstrated by comparing the temporal evolution of the total plasma content in 
the flux tube as seen from the two model. 

Figure 11 shows the total content as a function of time. As in Figure 10, for 
the hydrodynamic model the curves are for different values of the imposed velocity 
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at the boundary. It is seen that for V 0 = 0.39 V t , the hydrodynamic model yields 
nearly the same total plasma content as the kinetic model with nearly the same 
rate of increase in it for t < 12 hours. At later times, the content from the kinetic 
model shows a tendency toward saturation because the rate of increase in the content 
continuously decreases. Even in the hydrodynamic model there is a tendency towards 
the decreasing rate, but the decrease is much slower. This difference in the influx of 
the ions form the two model has a simple explanation. In the kinetic model, some 
of the ions have the liberty to exit the flux tube as they are scattered by Coulomb 
collisions, or as they simply flow out. On the other hand, in the hydrodynamic model 
the plasma entering the flux tube can leave the system only through the opposite 
boundary, where the flow velocity becomes exceedingly small after the shock phase 
(t > 2 hours). This implies that in the hydrodynamic model there is no provision 
for the plasma to leave the system. The slight tendency toward the saturation in 
the hydrodynamic model is due to the changing plasma condition near the boundary 
where the flow originates. As the plasma density near this boundary increases, the 
influx into the flux tends to decrease. 

Figures 10 and 11 also show the equatorial density and the total plasma content 
from the hydrodynamic model for V 0 = 0A\\ and V 0 = 0 . It is seen that even for 
Vo = 0 , the equatorial density and the total content are increasing with time and, 
in about 48 hours, they tend to approach the corresponding results from the kinetic 
model. It may sound strange how refilling can occur with a boundary condition of 
V 0 = 01 The refilling of a flux tube with zero flow velocity as boundary conditions 
in a hydrodynamic model was previously described by Singh et al.< [1986]. 

The ion flux developing near the boundaries from the kinetic model and the 
hydrodynamic model for V 0 = 0.39 Vi during the early times (t < 12 hours) is 1.8 x 10 
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ions cm~ 2 s~ l . In the kinetic model, this flux continuously decreases because some 
of the ions entering the flux tube eventually leave. However, in the hydrodynamic 
model the plasma entering the flux tube remains in it, causing the over- refilling for t 
> 12 hours when V' 0 = 0.39V*. When V 0 = 0. the hydrodynamic model yields a flux 
of < 10 8 ions cm- 2 s- 1 and it continuously decreases as the forces (determined by 
density and temperature gradients) on the ions accelerating them into the flux tube 
from the boundary cells diminishes with the refilling. It is worth mentioning that the 
comparison carried out above is based on the simplified boundary conditions in the 
kinetic and fluid treatments and a heuristic treatment of the heatflow. A comparison 
of the plasma treatments without these simplifications will be worthwhile. 

A comparison of plasma distributions in the flux tube at a relatively late time 
[t = 48 hours), as obtained from the two models, is shown in Figures 12a to 12d. For 
the hydrodynamic model, r/ = 0.3, and the distributions are given for three values of 
the boundary velocity, V 0 = 0.39 V„ 0.1 V t , and 0. Density profiles in Figure 12a 
show the over- refilling for V 0 = 0.39V t , but when V 0 is reduced below 0.lV t , the 
density profiles from the two models disagree near the boundary A = \ 0 , but away 
from it the agreement considerably improves. The disagreement near the boundary 
is also reflected in the velocity profiles in Figure 12b. Despite the above disagreement 
in the density and velocity profiles near the boundary, the temperature structures 
obtained from the two models are nearly identical. Temperature is nearly isotropic 
( T\\ = Ti); it rises to about 0.4 eV in the equatorial region from the boundary value 
of 0.3 eV. In the late stage of the refilling, the similarity between the temperature 
profiles, despite the differences in the density and velocity profiles from the two 
models, can be understood by examining the temperature equations (3) and (4) and 
the flow properties. Since the flow velocity is small, the velocity profile has little 
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effect on the temperature profiles. The maximum flow velocity in a localized region 
near the boundary is 2.5 km/s, compared to the thermal velocity of 5.5 km/s and 
ion-acoustic speed of 10 km/s. The density distribution affects the temperature 
distribution through the heatflow terms in equations (3) and (4). In the late stage of 
the flow when the gradients have smoothed out and the densities are relatively large, 
the density distribution also has insignificant effects on the heat flow. 
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4. Conclusion and Discussion 


We have carried out a comparison between semi-kinetic and hydrodynamic mod- 
els for plasma flow along closed magnetic field lines. The comparison has direct 
relevance to the problem of plasmaspheric refilling. It is found that the compari- 
son does not depend only on the plasma physics afforded by the models, but it also 
strongly depends on the boundary condition on the flow velocity. In a kinetic model, 
an appropriate boundary condition is to prescribe the velocity distributions of the 
inflowing ions to be half- Maxwellian for V > 0 at A = \ 0 and for V < 0 at A = — \ 0 . 
In the hydrodynamic model, this boundary condition corresponds to the drift veloc- 
ity Vq = ( kTo /2nm } 1/2 . A comparison of results from the two models with such 
boundary conditions revealed the following important features of the flows. 

1. When supersonic flows develops in response to a sudden depletion in a 
flux tube, the hydrodynamic and kinetic models yield distribution of den- 
sity, flow velocity and temperatures in generally good agreement. The 
temperature distributions in the region of supersonic flows are found to 
be remarkably similar, showing small effect of the heatflow. It is worth 
pointing out that Demars and Shunk [1991] compared the behavior of 
a highly supersonic plasma flow from a hydrodynamic model based on 
a more complete (16-moment) set of equations with that from a semi- 
kinetic model, demonstrating a good agreement. We have demonstrated 
here that, for a highly supersonic flow, even a much simpler set of hydro- 
dynamic equations are adequate. It is physically explained by the fact 
that the transport of heat in a supersonic flow is dominated by the large 
drift velocity and not by the heatflow process. Mathematically speaking, 
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it implies that the heatflow terms are negligibly small compared to the 
convective terms in the temperature equations. 

2. Both models show reflection of the supersonic flow when it penetrates 
deep into the opposite hemisphere. Since even a two-stream hydrody- 
namic model can not handle the counterstreaming for a given flow, the re- 
flection automatically leads to a shock formation [Rasmussen and Schunk , 
1988; Singh , 1991]. The shock first moves upward toward the equator 
and then downward to the ionospheric boundary. An examination of the 
plasma conditions for shock formation shows that the shock seen in the 
hydrodynamic model is an artifact of the model; the ion beams are found 
to be too warm to excite the ion-ion instability which can subsequently 
produce a shock. The semi-kinetic model shows the development of coun- 
terstreaming for the flow; the counterstreaming advances to the equator 
and downward to ionospheric boundary. It turns out that the transit time 
of the shock all the way to the ionospheric boundary and the time for the 
counterstreaming to spread to this boundary are nearly the same, about 
2 hours. In a previous paper, Singh [1991] reported the shock transit time 
to be about 4 hours, which is in error due to a normalization factor of 2. 
In view of the short transit time of the shock, the shock formation does 
not significantly affect the refilling as evidenced by the comparison of the 
flows from the hydrodynamic and kinetic models for later times. 

Lack of shock formation in the equatorial region, when the ion beams 
begin to interpenetrate [ Banks et al.< 1971] is uncertain in view of the 
spatial and temporal resolutions afforded by a large-scale model and the 
simplicity in handling the electron dynamics by the Boltzmann law. 
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3. After about 2 hours, the flow in each hemisphere becomes subsonic with 
respect to the ion-acoustic speed. This is seen from both the models. 

4. A comparison of the total plasma contents and the equatorial densities 
from the two models indicates a good agreement up to about t ~ 12 hours, 
after which the hydrodynamic model indicates over-refilling of the flux 
tube. The over-refilling is traced to the inability of our hydrodynamic 
model to control the net plasma inflow by the returning particles. The 
inflow is determined by the imposed boundary conditions and the outflow 
of plasma is exceedingly small. On the other hand, in the kinetic model 
the influx gradually decrease due to ions returning from the flux tube, 
showing a tendency toward saturation in the refilling in about 2 days. It 
is worth pointing out that it will be useful to perform a study comparing 
the models based on the kinetic and hydrodynamic treatments by relaxing 
some of the simplifications in terms of boundary conditions and in handling 
of the heat flow in the latter treatment. The boundary conditions can be 
relaxed by including the ionospheric plasma generation processes at low 
altitudes [ Gaiter and Gombosi , 1990]. 

When the boundary flow velocity in the hydrodynamic model reduces below 
V 0 = {kT 0 j 2irm ) 1/2 , there is an initial underfilling, but eventually the refilling from 
this model catches up to that given by the semi-kinetic model. For example when 
V 0 = 0, the degree of refilling from the two models, in terms of both the equatorial 
density and the total plasma content in the flux tube, becomes approximately the 
same in about 2 days. 

In some previous studies [Singh, et ai, 1986; Rasmussen and Schunk , 1988; Singh, 
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1991], a boundary condition of zero flow velocity was used. It may appear strange 
that a refilling occurs with this boundary condition on the flow velocity. The issue 
is briefly revisited here. 

From the comparison of the plasma contents and the equatorial densities given 
by the models, it is concluded that after about 12 hours, the choice of boundary 
condition in the hydrodynamic model is quite uncertain . In view of this uncertainty, 
the choice of zero-velocity boundary could be useful during the late stage of the 
refilling, it yields under-refilling only near the boundaries, where the density and 
average flow velocity show a discontinuity in the flow. Otherwise, over the rest of the 
flux tube the density and flow velocity are in quite good agreement with those given 
by the kinetic model. 

The hydrodynamic model described here is a two-stream model and includes equa- 
tion for the parallel and perpendicular temperatures. Single-stream hydrodynamic 
models [Singh et al., 1986; Guiter and Gombosi , 1990] suffer from the shortcoming 
that they generate shocks at the equator whether the plasma conditions allow them 
or not. The single- and two-stream models with assumed temperature isotropy suffer 
from the shortcoming that the shock transit time is fairly long and major part of 
the refilling occurs through supersonic flows from the ionospheres [ Singh et al., 1986, 
Rasmussen and Schunlc. 1988; Singh, 1991]. This is in contrast to the two-stream 
model with a self-consistent treatment of the temperature anisotropy; this model 
yields evolution from supersonic to subsonic flows at the same time scale as the ki- 
netic model. In this sense, the heuristic treatment of the heatflow described in this 
paper appears to be adequate. This treatment also appears to be adequate even in 
the subsonic stage as long as the flow velocity near the boundary is relatively large 
near the thermal speed, for example, for t < 12 hours in our present calculations. 
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However, when the flow velocity becomes sufficiently low so that a large fraction of 
injected ions in the kinetic model begin to return from the immediate vicinity of the 
boundary, the boundary conditions for the hydrodynamic model diverges from that of 
the kinetic treatment, because this model does not allow for a return flux for a given 
stream. Can this situation be improved by a more rigorous treatment of the heatflow 
and/or by properly including the ionospheric plasma supply [ Guiter and Gombosi , 
1990]? In order to answer this question, it will be useful to compare models based on 
(1) the heuristic heatflow treatment, (2) a more sophisticated treatment of the heat- 
flow using 16-moment set of transport equations, and (3) the semikinetic treatment, 
and all models properly including the analogous ionospheric boundary conditions. If 
the model (1) compares well with the latter ones, the computational effort in using 
transport equations for modeling space plasma will be considerably reduced. 
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Figure Captions: 


Figure 1 . Geometry of a closed magnetic flux tube. The latitudinal angle A and the geocentric 
distance r are shown. The ionospheric boundaries are at 3 = 0 (A = A 0 ), and s = 

‘^max(A — A o ) 

Figure 2 . Temporal evolution of the flow originating from the northern hemisphere; phase-space 
{s _ VJ|) plots are shown. The plot at t = .003 hour nearly shows the initial plasma 
in the flux tube. 

Figure 3 . Comparison of flows from semi-kinetic and hydrodynamic models at t = 30 minutes 
For the latter model flows are given for 77 = 0, 0.05 and 0.3: (a) density, (b) flow 
velocity, (c) parallel temperature and (d) perpendicular temperature distributions. 
For most latitudes, the curves are so close together that it is difficult to distinguish 
them. 

Figure 4 . Same as Figure 3 , but at t = 1 hour The hydrodynamic curves are distinguished by 
the presence of a shock which is manifested by sudden jumps in density, flow velocity 
and parallel temperature near A ^ —25. Shocks are indicated by the arrows. 

Figure 5 . Propagation of the shock is shown through the jump in the density profiles at (a) 
t hour, (b) f = 1.5 hours and (c) t — 2 hours. For the purpose of comparison 
the curve from the kinetic model and three curves for 77 = 0,0.05, and 0.3 from the 
hydrodynamic model are shown. Shocks are indicated by arrows. 

Figure 6 . Distribution of flow parameters from the kinetic model. The average flow velocity 
V 5 , ion-acoustic speed C s and the temperature ratio T\\/T e for ions with X j| > 0 are 
shown for the purpose of instability analysis at t = 30 minutes. The corresponding 
curves for Vjj < 0 can be deduced from the symmetry considerations. 
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Figure 7. Same as Figure 6. but at t — 1 hour. 

Figure 8. Same as Figure 4. but at t — 4 hours. 

Figure 9. Same as Figure 4, but at t = 12 hours. 

Figure 10. Temporal evaluation of the equatorial density from the kinetic model and from the 
hydrodynamic model for three values of Fo, 0.39 F t , 0.1 F f , and OV^. 

Figure 11. Temporal evolution of the total plasma content in the flux tube for the kinetic and 
hydrodynamic models. 

Figure 12. Distribution of (a)density, (b) flow velocity, (c) parallel temperatures and (d) per- 
pendicular temperature at t = 48 hours . All vertical scales are linear in this Figure. 
There are three curves from the hydrodynamic model for Vg = 0, 0.1 1 < . and 0.39 V t . 
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